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PACS 64 . 70 . D- - Solid-liquid transitions 

PACS 61 . 72 . J — Point defects and defect clusters 

PACS 83. 10. Rs - Computer simulation of molecular and particle dynamics 

Abstract. - Using molecular dynamics simulation, we study structural and dynamical hetero- 
geneities at melting in two-dimensional one-component systems with 36000 particles. Between 
crystal and liquid we find intermediate hexatic states, where the density fluctuations are en- 
hanced at small wave number k as well as those of the six- fold orient ational order parameter. 
Their structure factors both grow up to the smallest wave number equal to the inverse system 
length. The intermediate scattering function of the density S(k, t) is found to relax exponentially 
with decay rate Tk oc k z with z ~ 2.6 at small k in the hexatic phase. 



0$ Introduction. — Since a simulation by Alder and 
£h Wainwright [1], much attention has been paid to the two- 
i dimensional (2D) melting in simple one component parti- 
^3cle systems [2]. However, it has been controversial whether 
gthe transition is first order as in three-dimensional melt- 
Qing [1-5] or is continuous as predicted by Halperin and 
•—Nelson [6]. They presented a defect-mediated melting 
mechanism and a "hexatic phase" in a temperature (or 
^density) window between crystal and liquid. In the hex- 
(NT)atic phase, the bond-orientation correlation function ge(r) 
VOof a sixfold orientation order parameter ipe(r) decays alge- 
braically, indicating a quasi- long- range orientational cor- 
relation. Afterwards their prediction has been confirmed 
Oin experiments [7-12] and in simulations [13-20]. In the 
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hexatic phase, defects have been observed to proliferate 



^^with increasing the temperature T or decreasing the den- 
• *sity n. As other theories, Chui proposed a melting mech- 
.^anism mediated by grain boundaries [21], while Saito ar- 
^ gued that the 2D melting can be either continuous or first 
Unorder depending on the specific details of the system [15]. 
As a marked feature, a number of experiments and sim- 
ulations [1-3,7-14,22] have observed heterogeneities in the 
hexagonal structures and in the particle trajectories, de- 
veloping around the transition. In particular, appreciable 
dependence on the system size has been encountered in 
the calculations of the equation of state [1,5,17,20] and 
the local fluctuations of ^(r) [5,20]. However, the het- 
erogeneities in the 2D melting have not yet been well un- 
derstood. In this Letter, we will visualize them using a 



disorder variable representing deviations from the hexag- 
onal order [23] and bond breakage used in analyzing glass 
dynamics [24]. It is also a fundamental issue whether the 
isothermal compressibility Kt = (dn/dp)T /n remains fi- 
nite or tends to infinity in the hexatic phase. In simula- 
tions [17,20], the pressure p was a weakly decreasing func- 
tion of n, apparently suggesting (dp/dn)T < 0, in the hex- 
atic phase. Hence, we will calculate the structure factor of 
the density S(k) at small wave number k to see whether 
the thermodynamic relation lim^^o S(k) = n 2 TK T holds 
or not. The intermediate scattering function S(fc, t) will 
then emerge as a new informative quantity. 

Numerical Method. — Our 2D system is com- 
posed of N = 36000 particles interacting via a truncated 
Lennard- Jones potential of the form 



v(r)=4e [(a/r) 12 - (a/r) 6 ] - C, 



(1) 



which is characterized by the energy e and the range a. 
For r > r cut = 3.2<j, we set v(r) = with the constant 
C ensuring the continuity of v(r) at the cut-off. The sys- 
tem volume V is kept fixed such that (j) = Na 2 /V = 0.9. 
Then the system length is given by L = V 1 / 2 = 200<r. 
We integrated the equations of motion using the Stormer- 
Verlet algorithm (a sort of the leap-frog method) under 
the periodic boundary conditions using the Nose-Hoover 
thermostat. The time step of integration is 0.002r with 



(2) 
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Fig. 1: Bond-orientational correlation function g^ir) on a log- 
arithmic scale for T = 1.0, 1.1, 1.15, 1.2, 1.25, 1.3, and 1.4 from 
above, obtained from data averaged over long-time and over 
runs starting with 5 independent initial conditions. Below the 
curves of T — 1.1, 1.15, and 1.2, dotted line segments are guides 
to eye with slopes —0.03, —0.09, and —0.36, respectively. 

m being the particle mass. In our simulations, we first 
quenched the system from T = 2e/kB to 0.2e//c# into a 
crystal state. After a relaxation time of 5 x 10 3 r, there was 
no appreciable time evolution in physical quantities such 
as the pressure and g^ir). We then increased T to a desired 
value. The time t is set equal to at this temperature 
increase. We continued the simulation until t = 2 x 10 4 r. 
Hereafter we will measure space, time, and T in units of 
(j, t, and e/ks, respectively. 

In 2D dense particle systems, a large fraction of the 
particles are enclosed by six particles and the local order 
is represented by the sixfold orientation [6]. We define 
the orientation angle aj in the range [— 7r/6,7r/6] for each 
particle j at position rj using the complex number 

Yl exp[6^ fc ] = |* i |e 6 ^, (3) 

fcGbonded 

where the summation is over particles "bonded" to the 
particle j. The two particles j and k are bonded if \rj — 
rk\ < 1.25a. Ojk is the angle between Vj — and the 
x axis [23,24]. Next we construct another nonnegative- 
definite variable representing the degree of disorder for 
each particle j by [23] 

D j= 2 [1- cos 6(^-^)1. (4) 

/cGbonded 

Here Dj is nearly zero for a perfect crystal, but is large in 
the range 5 — 20 for particles around defects. Thus Dj is 
convenient in visualizing the structural inhomogeneity. 

In terms of aj the sixfold orientation order parameter 
iPq(v) is defined as [6] 

N 
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(c)T=1.2 (d)T=1.3 
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Fig. 2: Snapshots of disorder variable Dj of 36000 particles at 
t = 1.2 x 10 4 for (a) T — 1.0 (crystal), (b) T — 1.1 (hexatic), 
(c) T = 1.2 (hexatic), and (d) T — 1.3 (liquid). Colors are 
given according to the bar on the top. Particles with Dj > 4 
are written in red. 

In the hexatic phase, the bond-orientational correlation 
function decays algebraically as 

g e (r) = (Mr)M0y)~r- r >, (6) 

where r = |r|. Theoretically [6], the exponent r] depends 
on T and n = N/V in the range < r] < 1/4. 

Structural Heterogeneity. — In fig. 1, the curves 
of ge(r) are displayed, which are the long-time averages 
over the snapshots produced by 5 independent runs in the 
range 10 4 <t<2x 10 4 . In the hexatic phase at T = 1.1, 
1.15, and 1.2, the exponent r] in eq. (6) is 0.03,0.09, and 
0.36, respectively, where 0.03 at T = 1.1 is very small and 
0.36 at T = 1.2 even exceeds the theoretical upper bound 
1/4. The hexatic-liquid and crystal-hexatic transitions are 
continuous, taking place at T = 1.2 and 1.0 respectively. 
For our system size, it is still difficult to determine the 
transition temperatures precisely. For T > 1.25 the sys- 
tem is in liquid with g&(r) decaying exponentially, while 
for T < 1.0 the system is in crystal without appreciable 
decay of g 6 (r). 

In fig. 2, we display snapshots of Dj of all the parti- 
cles in a crystal phase at T = 1.0 (a), in hexatic phases 
at T = 1.1 (b) and 1.2 (c), and in a liquid phase at 
T = 1.3 (d). In fig. 3 (a), a more expanded snapshot 
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Fig. 3: Snapshot of disorder variable Dj (a), local packing density pj (b), and particles with neighbor number different from 
six (c) for 36000 particles at T — 1.15. These three panels exhibit heterogeneous patterns stemming from the same long-range 
structural disorder. 
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Fig. 4: Structure factors (a) Se(k) and (b) S(k) for the six- 
fold orientation order and the density, respectively, for T = 
0.8,1.0,1.1,1.15,1.2, and 1.4, calculated as the long-time av- 
erage. Averages are taken also over 5 independent runs. Use 
of mark in (b) is similar to (a). They are enhanced at long 
wavelengths in the hexatic phase. 



of Dj is given at T — 1.15. The average disorder param- 
eter D = Y,j D j/ N is im , L 72 > 2 - 12 > 2 - 55 > and 332 for 
T = 1.0, 1.1, 1.15, 1.2 and 1.3 respectively. In the hexatic 
phases (1.1 < T < 1.2), marked heterogeneity emerges 
on large scales among crystalline and disordered regions, 
though there are no sharp boundaries. The patterns are 
fractal-like, resembling the critical fluctuations near the 
Ising criticality. Also shown in fig. 3 are (b) the local areal 
density pj (to be defined below) and (c) the particles with 
neighbor number being different from six. A common par- 
ticle configuration was used for these three panels. Here 
vj = pj 1 is the volume of particle j in the Voronoi cell con- 
struction. We treat its inverse pj as the local density at the 
position of particle j. Its variance V (= (J2j(Pj ~ P) 2 /N) 
is about 0.05 with p = . Pj/N = 0.9 here. Comparing 
the two panels we can see that the particles with larger Dj 
tend to have smaller pj. In the literature [7,10,12], defects 
have been detected around the particles with five or seven 
neighbors, so we used this Voronoi method to produce 
fig. 3 (c). It exhibits essentially the same heterogeneity as 
that of Dj (a), though only discrete particles are selected. 
However, near melting, using Dj is more quantitative and 
appropriate to characterize the diffuse disordered regions 
extending on large scales. 

The algebraic decay of g^ir) in eq.(4) arises from the 
heterogeneity in figs. 2 and 3. Furthermore, fig. 3 (b) 
indicates that small density differences exist among the 
crystalline and disordered regions. In fig. 4, we thus show 
the structure factors of the hexagonal order and the den- 
sity 

S 6 (k) = J dre ik - r g 6 (r) = (\^ 6k \ 2 ), (7) 
S(k) = Jdr e ikr (Sn(r)Sn(O)) = (|n fc | 2 ), (8) 

where Sn(r) = ^ S(r — Vj) — n is the microscopic density 
deviation. The V^fc and are the Fourier components of 
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Fig. 5: Pressure p (right axis) vs density n at T = 1.15 with 
a small negative slope in the hexatic phase. Compared are re- 
sultant nT/(dp/dn)T and S(k m in) (left axis), which coincide 
outside the hexatic region confirming the the compressibility 
relation. In the hexatic phase, (dp/dn)T is apparently nega- 
tive and the long wavelength limit S(0) is not attained for our 
system size. 

ijjQ (r) and n(r). These structure factors are the averages 
over the angle of the wave vector k = (k Xl k y ). The small- 
est wave number /c m i n is defined as tt(1 + V2)/L = 0.035, 
where L = 200 is the system length. The structure 
factors at k = k m { n are the averages of the data at 
k = 2ttL- 1 (±1,±1) and 2ttL- 1 (±v / 2, ±\/2). In fig. 4(a), 
the growth Sq (k) = AQk~ 2+r] can be seen at small k in the 
hexatic phase at T = 1.1, 1.15, and 1.2 in accord with eq. 
(6), where A 6 = 0.53 and r] = 0.09 at T = 1.15. In fig. 4 
(b), S(k) grows at small k in the hexatic phase, but its 
amplitude is very small and S(k m { n ) remains smaller than 
the peak height at k ~ 2tt by two orders of magnitude (see 
the inset). In fact, at T = 1.15, its curve may be fitted 
to S(k) = 0.017 + 1.23 x lO- 3 /^ 1 - 55 for k < 1. The small 
coefficient (~ 10 -3 ) here arises from small density differ- 
ences among the crystalline and disordered regions. For 
our system size, these structure factors do not saturate 
even at k = k m [ n in the hexatic phase. 

If S(k) saturates to a long wavelength limit 5(0) = 
lini/^o S(k) in the thermodynamic limit L — > oo at 
fixed density, the compressibility is given by Kt = 
(dn/ 8p)t /n = S(0)/n 2 T. From our simulation only, how- 
ever, we cannot exclude the possibility of S(k) — > oo (as 
k — > 0) in the hexatic phase, where (dp/dn)T = ulti- 
mately holds in the thermodynamic limit. As in fig. 5, we 
also performed simulations with 36000 particles by vary- 
ing the volume V, where shown is the pressure p (av- 
erage of its microscopic expression) vs the density n at 
T = 1.15. Outside the hexatic region, the long wave- 
length limit S(0) is attained and the compressibility rela- 
tion S(0) = nT/(dp/dn)T surely holds. Here, as in previ- 
ous work [20] , p apparently exhibits a small negative slope 
in the hexatic density range. Thus we need to use much 
larger system sizes to settle this issue. In such simulations, 
the long wavelength fluctuations (with k < 10 -2 ) need to 
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Fig. 6: (a) Mean-square displacement vs t for various T on a 
logarithmic scale, yielding the diffusion constant D s . In crystal 
at T = 0.6 and 0.8, its linear increase after a plateau is due 
to defect motion, (b) Surviving-bond number Nb(t) vs t on a 
semi-logarithmic scale, decaying as e _t ^ Tb . Here D s ~ 1/t&. 

be equilibrated on extremely long times. 

Dynamics. — The dynamics has not yet been well 
studied at the 2D melting. For the particle displacement 
Arj(t) = Vj(t + to) — Vj(to) in time interval £, fig. 6 
(a) displays the mean square displacement, ([Ar(t)] 2 ) = 
^Zj([Arj(t)] 2 ) /A/", which is the average over all the par- 
ticles and over the simulation time. The linear behavior 
([Ar(t)] 2 ) ^ 4D s t can be seen at long times for T > 0.6. 
The diffusion constant D s thus obtained increases as 0.181, 
1.13, 3.32, 8.23, 60.5, and 188 for T = 0.6,0.8,0.9,1,1.1, 
and 1.2, respectively, in units of 10" -v-v- 1 . In crystal 
(0.6 < T < 1.0), the curves exhibit a plateau followed by 
the linear growth. Similar two-step behavior is well-known 
in supercooled liquids [24], but it is here due to motions 
of defect clusters composed of several particles with finite 
Dj (see fig. 2 (a) ). Such clusters were observed in 2D 
colloidal systems [12]. In the hexatic phase, on the other 
hand, the plateau disappears and D s grows abruptly. 

To examine the particle-configuration changes, we in- 
troduce the bond breakage [24]. For each particle config- 
uration at a time to(~ 10 4 ) after long annealing, a pair of 
particles i and j is considered to be bonded if 

ri j (t ) = \r i {to)-r j {t )\<A 1 , (9) 
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Fig. 7: Snapshots of broken bonds in three consecutive time intervals [to + (£ — 1) At, to + £At] {£ = 1,2,3) for (a) At = 800 at 
T = 0.8, (b)At = 100 at T = 1.0, and (c) At = 10 at T = 1.2. Colors: ultramarine (£ = 1), magenta (£ = 2), and black (I = 3) 
in the chronological order. Corresponding snapshots of Dj of (a) and (b) at £ = 3 are shown in fig. 2 (b) and (c) 



where we set A\ = 1.2 (around the peak distance of the 
pair correlation function). After a time interval t, the 
bond is regarded to be broken if 

ryito+t) > A 2 , (10) 

where we set A 2 — 1.6. In fig. 6 (b), we plot the number 
of the surviving bonds N b (t) vs t for various T. It is equal 
to the initial bond number N b (0) (~ 8.4 x 10 4 ) minus 
the number of the broken bonds. It may fairly be fitted 
to the exponential form e~ l l Th . The bond life time r b is 
determined by N b (r b ) = N b (0)/e. Then r b = 91, 8.1, 
1, 0.3, and 0.09 for T = 0.6,0.8,1, and 1.2, respectively, 
in units of 10 3 r. Here we notice that the product D s r b 
is between 1 and 2, where D s is the diffusion constant 
determined in fig. 6 (a). Thus 

D.-Tt 1 , (11) 

which demonstrates that the particle motions are caused 
by the configuration changes in our jammed states. 

We examine how the dynamic heterogeneity evolves in 
time. In fig. 7, we pick up the particles with broken bonds 
in three consecutive time intervals and mark them in ul- 
tramarine, magenta, and black in this order. Here we set 
(a) At = 800 at T = 0.8, (b)At = 100 at T = 1.0, and 
(c)At = 10 at T = 1.2. In crystal (a), the evolution is 
due to the motions of defect clusters taking place in the 
form of string-like trajectories in each event. Such tra- 
jectories accumulate to form large-scale dynamic hetero- 
geneity on long time scales, though the defect number is 
small at each time in crystal. This picture explains the 
two-step behavior of ([Ar(t)] 2 ) in crystal in fig. 6 (a). In 
addition, two straight lines arising from dislocation glid- 
ing can be seen in (a), but such slip motions are rare in 
our system. In (b), the system is still in crystal, but the 
time scale of bond breakage is much faster. In the hexatic 
phase (c), the particles in the disordered regions are rel- 
atively mobile, while those in the crystalline regions are 



nearly immobile. The mobility of the particles is strongly 
correlated to the disorder variable Dj visualized in fig. 2 
(c). The dynamic heterogeneity can be seen over a rather 
broad temperature range around the melting, where the 
time scale changes dramatically [9,12-14,22,23]. 

Since the large-scale density fluctuations are enhanced 
as in fig. 4 (b), we are interested in their relaxation. As 
shown in fig. 8, we calculated the total intermediate scat- 
tering function 

S(M) = J dr e ik r (5h(r,t)5n(O,0)) (12) 

in the hexatic phase at T = 1.15 where the initial value 
S(k) is enhanced as in fig. 4 (b). In an early stage it 
undergoes an oscillatory decay arising from the acoustic 
propagation. We find that it then decays exponentially as 

S(M) = S(k)A k e~ Tk \ (13) 

for small « 1. The amplitude Ak approaches unity 
and the decay rate behaves as ^ F with z ~ 2.6 for 
small k. The dotted lines in fig. 8 represent this expo- 
nential form, which are excellently fitted to the numerical 
data. This slow decay arises from the evolution of large- 
scale density heterogeneities produced by the structural 
fluctuations. In the hexatic phase, the diffusion constant 
depends on the wave number as D k oc k z ~ 2 if it is intro- 
duced by Dk = Tk/k 2 . This is analogous to the thermal 
diffusion constant at the gas-liquid criticality, where z = 3 
in three dimensions [26]. However, we cannot explain this 
exponential relaxation in the hexatic phase at present. In 
passing, we also calculated S(k, t) in crystal and liquid 
outside the hexatic temperature window, where the den- 
sity fluctuations are much suppressed at small k. There, 
the long wavelength relaxation of S(k,i) is due to hydro- 
dynamic thermal diffusion with 1^ oc k 2 for small k. 
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Fig. 8: Normalized intermediate scattering function 
S(k,t)/S(k) in the hexatic phase at T = 1.15 for vari- 
ous small k. After the initial acoustic damped oscillation, it 
can be fitted to the exponential decay of the form Ake~ Fkt 
(solid lines) with Tk oc k 2 ' 6 . 

Summary and Comments. — We have examined 
the heterogeneities in structure and dynamics at two- 
dimensional melting of one-component systems at fixed 
volume V. In terms of the disorder variable Dj, struc- 
turally heterogeneous patterns have emerged unambigu- 
ously among the crystalline and disordered regions in the 
hexatic phase in figs. 2 and 3. Though very weak, we 
have noticed the presence of large-scale density fluctua- 
tions. As a result, at small wave number /c, the structure 
factor Sq (k) of the sixfold orientation order grows strongly 
as in fig. 4 (a), while the structure factor S(k) of the den- 
sity grows similarly but its amplitude is very small as in 
fig. 4 (b). Dynamically heterogeneous patterns have been 
obtained in fig. 7 in the crystal and hexatic phases. Crystal 
states are dynamically heterogeneous on long time scales 
due to relatively rapid motions of defect clusters. In the 
hexatic phase, dynamics is heterogeneous on short time 
scales, where the particles with high Dj tend to be mobile 
than those with small Dj. We have also calculated the 
intermediate scattering function S{k 1 t) 1 which has turned 
out to relax exponentially with the dynamic exponent z 
about 2.6 at small k in the hexatic phase. 

We make further comments, (i) In our simulation, meso- 
scopic coexistence of the ordered and disordered regions is 
realized dynamically in the hexatic phase. There are no 
sharp boundaries between the two regions and the free en- 
ergy penalty due to the structural inhomogeneity should 
be very small. These structural fluctuations resemble the 
critical fluctuations in Ising systems. It remains puzzling 
whether or not the long wavelength limit of S(k) tends to a 
finite constant in the thermodynamic limit. If so, it follows 
a finite compressibility, (ii) Growth (shrinkage) of the dis- 
ordered regions gives rise to an increase (a decrease) in the 
pressure p at fixed volume. We should perform constant- 
pressure simulations also to examine whether or not the 
hexatic phase exists depends on the boundary condition. 



(iii) In binary mixtures the size dispersity serves to pin the 
particle motions and the relaxation times are much longer 
than in one-component systems, (iv) Shear flow has been 
applied to glassy and polycrystal systems [24-26]. It is 
intriguing how applying shear can affect the hexatic state 
and the scenario of the 2D melting. 

* * * 
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